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The aim of this mini review is to survey the literature on the the study of the nonequilibrium 
dynamics of Fermi superfiuids in the BCS and BEC limits, both in the single channel and dual 
channel cases. The focus is on mean field approaches to the dynamics, with specific attention drawn 

■ to the dynamics of the Ginzburg-Landau order parameters of the Fermi and composite Bose fields, 
as well as on the microscopic dynamics of the quantum degrees of freedom. The two approaches are 
valid approximations in two different time scales of the ensuing dynamics. The system is presumed 

£N| ■ to evolve during and/or after a quantum quench in the parameter space. The quench can either 

be an impulse quench with virtually instantaneous variation, or a periodic variation between two 
values. The literature for the order parameter dynamics, described by the time-dependent Ginzburg- 
Landau equations, is reviewed, and the works of the author in this area highlighted. The mixed 

■ phase regime in the dual channel case is also considered, and the dual order parameter dynamics 
of Fermi-Bose mixtures is also reviewed. Finally, the nonequilibrium dynamics of the microscopic 

, degrees of freedom for the superfluid is reviewed for the self-consistent and non self-consistent cases. 

The dynamics of the former can be described by the Bogoliubov de-Gennes equations with the 
equilibrium BCS gap equation continued in time and self -consistently coupled to the BdG dynamics. 
The latter is a reduced BCS problem and can be mapped onto the dynamics of Ising and Kitaev 
models. This article reviews the dynamics of both impulse quenches in the Feshbach detuning, as 
well as periodic quenches in the chemical potential, and highlights the author's contributions in this 
F area of research. 
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I. INTRODUCTION 



Studies of the coherent quantum dynamics of closed many body systems have gained considerable interest in the 
last two decades. This interest has been brought about by such states becoming experimentally accessible using 
ultracold atoms, cooled to near absolute zero temperatures by the use of optical molasses [l[ , evaporative cooling 
laser culling [3fl , cavity QED methods [4| , and other experimental techniques. The experimental realization of the 
Bose Einstein condensate by Cornell, Ketterle and Weiman in 1995 @ , as well as the realization of the Bosonic 
Superfluid Mott-Insulator transition by Greiner et. al. in 2002 [5] (theoretically predicted by Jaksch et. al. in 
1998 [6] ), motivated strong theoretical interest in the coherent dynamics of ultracold bosons. Previously, the regimes 
, where fermionic coherent dynamics can be realized were unattainable in either solid state systems or with ultracold 
atoms. The difficulties for ultracold atoms were eventually circumvented by experiments performed over the first 
\Q \ decade of the 21 st century. Studies of fermionic atoms gained momentum after these experiments realized pure 
fermionic BCS states Q , as well as the BCS-BEC crossover that was predicted theoretically a decade earlier Q . 
ON ' The formation of Fermi superfiuids were realized by the use of strong four fermion interactions caused by tuning a 
. homogeneous magnetic field close to the Feshbach resonance of the confined atoms @ . Such strong interactions are 
t— I 1 difficult to generate with bosonic atoms due to three body losses. However, three body losses are significantly lower 
with fermions due to the Pauli principle. Theoretical interest in Fermi superfluid dynamics paralleled experiments 
involving dynamical phenomena such as trap expansions, collective excitations, and vortex dynamics. Theoretical 
studies of the exact dynamics focused on BCS- BEC systems near unitarity or in the deep BCS regime. However, the 
■ primary theoretical interest continues to remain with Bosonic dynamics, particularly the Bose-Hubbard model, and 
Ising-like spin models. 

This mini review focuses on theoretical works investigating the other end of the spectrum i. e, the dynamics of Fermi 
superfiuids. The exclusive focus is on mean field dynamics at T — as the system is quenched in various ways. The 
literature can be broadly divided into two kinds, one that looks at the macroscopic dynamics of the order parameters, 
covered in section [Til and one that looks at the dynamics of the microscopic degrees of freedom, covered in section Hill 
In the former case, the excited quasiparticles in the Bogoliubov de-Gennes spectrum of the superfluid relax very 
quickly compared to the relaxation of the time dependent superfluid gap A(i), and excitations of the quasimolecular 
states relax similarly [lOL |12[ . Thus, the dynamics of the system can be described by a mean field treatment of the 
macroscopic order parameter A(t) in the single channel Feshbach resonance, and a two-parameter system spanned 
by A(i) and the BEC matter field bo(t) in the dual channel case. In the latter case, the Bogoliubov excitations 
relax much more slowly, and the system must be described by the dynamics of the quasipa rticles with a coupling 
to a self-consistent equation for A(t). These regimes are discussed in more detail in |13l - fl7| . as well as in [l8l. fl9|. 
This review restricts attention to the mean field dynamics of population balanced Fermi gases at T = for Fermi 
superfiuids that can be modeled by integrable closed quantum systems. The dynamics of fluctuations, nonintegrable 
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systems, and systems connected to thermal reservoirs, involve diagrammatic techniques (see, for instance, ref [161 ]) 
, time dependent Density Matrix Renormalization Group techniques (see, for instance, |20] ), and the Eigenstate 
Thermalization Hypothesis (see, for instance, ref respectively. These are beyond the scope of this review, and 

the reader is directed to the references cited above for details. 

This mini-review article contains results from [Hill HHU. Section OH contains material from [24j that was 
presented by the author at the 2012 National Conference on Nonlinear Systems and Dynamics (Quantum Chaos 
minisymposium), held at the Indian Institute of Science Education and Research (USER), Pune, India. Its new 
contribution consists in generalizing Il4| to explicitly time dependent quenches, as well as suggesting future directions 
in this area of research based on [l4l.ll7l.l24j . This mini review also connects the different models and paradigms used 
in the mean field dynamics of Fermi BCS systems, and highlights literature in an area of research that has gained 
considerable interest in the last few years. 

II. DYNAMICS OF THE ORDER PARAMETERS 

This section details studies of the dynamics of the macroscopic order parameters that arise from the microscopic 
BCS-BEC Hamiltonian. Studies of these regimes have been motivated by the BCS-BEC crossover, where the renor- 
malized interaction between the Cooper pairs changes sign, and the size of the Cooper pairs decrease. These pairs 
eventually become diatomic molecules that obey Bose statistics and can form a Bose Einstein Condensate. The equi- 
librium equation of state of the macroscopic order parameters has also shown the BCS-BEC crossover [ll[ , and this 
section details works that involve the nonequilibrium extension of these equations. 

The dynamics of ultracold Fermi superfluids can be understood by approximate models based on how the time 
scales of the dynamics compare to the relaxation of the quasiparticles of the equilibrium BCS Hamiltonian. The 
order parameter relaxes at a rate given by ftA -1 , where A is the macroscopic BCS energy gap responsible for 
superconducting/superfluid behavior. When this relaxes much more slowly than the energies of the quasiparticles, 
they reach their equilibrium quickly and the nonequilibrium state can be characterized by the macroscopic order 
parameters. This, however, only occurs in very specific regions in the full parameter space [Til ] near the critical 
temperature. More general mean held treatments are discussed in section 1X111 . Nonetheless, recent works have focused 
on investigating the order parameter dynamics of Fermi Bose mixtures. The most general dual channel model has been 
espoused by Holland et. al. [25[ and Timmermans et. al [26j as a model for BCS superfluidity of neutral fcrmionic 
atoms that interact attractively via a Feshbach resonance that is deep enough to induce molecular bonding. Most 
recently, the ensuing dynamics of mean held order parameters following an impulse quench of a system parameter at 
t = 0— has also been studied. (Roy [22|). 

The dynamical equations of motion can be obtained from the D— dimensional Timmermans' Hamiltonian (2|| H tm 
for a Fermi-Bose mixture at T = confined in a unit volume. Large impulse quenches in the Feshbach detuning can 
produce persistent oscillations in the matter wave field due to the nonlinear modes, similar to the collapse and revival 
of matter fields in such systems reported in the literature [ill [27[ . The Timmermans ' Hamiltonian in momentum 
space can be expressed in terms of ak. CT , the annihilation operator for the BCS fermions for momentum k and (pseudo) 
spin a =t, |, and 6 q , the operators of the composite bosons that are formed by the Feshbach attraction of two fermions 
of momenta ±p + q/2 and opposite spin. This yields [22| 

H tm = X ( £k ~ Mf) aL a k CT - I'M zl a kt a -k| a k';«k't 

k.tr kk' 

+ X (K + 2V - Pb) kqfrq + U B Y. & q' 1 5 qi 6 q2^qi^qi+q 2 ,qi+qi 

q qiqiqaqi 

+ 9r X ( 6 q a P+§t a -p+§ I + h - C ) ■ C 1 ) 
kq 

Here, the first line represents the fermions in momentum space with the BCS four-fermion interaction. The free fermion 
energies are represented by e^. The second line represents the Hamiltonian for the closed channel quasimolecular 
bosons 6 q (& q ) with energies E^ + 2v — hb and point contact repulsion of amplitude ub- Finally, the last line represents 
the atom-molecule coupling, through which the Fermi and Bose condensates exchange momentum q. Huang et. al. [l2j 
have obtained the excitation spectra of this Hamiltonian in the mean field, and Roy [22j has demonstrated that the 
relaxation time of these excitations are faster than E~ x , where E c — 2 \/ (i 2 ? + | A + g r bo\ 2 , A is the equilibrium BCS 
gap computed from the gap equation [28|, and bo is the ground state population of the bosons. Therefore, impulse 
quenches that are slower than E^ 1 allow an analysis of the dynamical system based only on a ground state BEC 
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field and fermion pair field to completely capture the non-equilibrium state. The Timmermans Hamiltonian can be 
rewritten in coordinate space as follows. 



H tm = / d D xxH tm (x), (2) 



where 



U tm {x) = ^2 yt>t( x ) i h F(r) - Hf] ^0*0 j - \uF\4>\(x)(l)l(x)(l) l (x)(f> 1 .(x) 



+ tf(x) [h B (r) + 2u~ fig] b(x) + u B b\x)b(x) [bHx)b(x) - i] 

+ g r [bHx)^{x)4> i {x) + h.c] , (3) 

the interactions of the excited quasimolecules have been ignored, and x = (v,t'). The fermion and composite boson 
fields in this system are represented by the operators <j> a {x) and b(x) respectively. Equation [3] describes a system of 
ultracold electrically neutral two-component fermions interacting attractively. The first line in equation [3] represents 

the Fermi-BC'S part of the Hamiltonian. Here, /ip(r)(/iB(r)) = — 2 ^) m + ^( r ) are the single particle Hamiltonians 

for the fermions (Bosons), and m is the fermion mass (the mass of the composite bosons is thus 2m). The second 
line represents the Hamiltonian of the composite bosons [HI, [26[ . Here, the Feshbach threshold energy (also called 
the Feshbach 'detuning' from the molecular channel to the continuum [26]) is represented by 2v, and ub represents 
the amplitude of the repulsion between the composite bosons. The final line describes the Feshbach resonance that 
leads to the fermions binding to (or dissociating from) the composite bosons, with g r representing the atom-molecule 
coupling. Note that the chemical potentials satisfy [is =2fip. 

The mean field dynamics of the Fermi and Bose order parameters can be obtained from the Timmermans Hamilto- 
nian via Ginzburg Landau Abrikosov Gor'kov theory, where the path integral grand partition function is transformed 
via the Hubbard Stratonovich method [28l [29j , introducing the fermion gap parameter A(t) via a Gaussian identity, 



J T>[A*, A]exp - J dx 



,A*(x')A(x'Y 

- 1. (4) 



u F 



Here, T>[A*, A] is the path integral measure of the gap field. Performing the Hubbard-Stratonovich transformation, 
A — > A + |uf|<^|,(^|- and A* — > A* + \up\(f>^(j)^, cancels out the four-fermion term from H tm . Now, integrating out the 
fermion fields using formal Grassman calculus allows for an effective action in terms of fermionic and bosonic order 
parameters. Machida and Koyama (30j obtained the dynamical equations of the coupled Fermi and Bose fields using 
the method summarized above for the special case of free particles, where the excitations of order E c were taken into 
consideration. This yielded spatio-temporally varying fields given by the solutions to the following dynamical system. 

i^-A(x) = aA(x)+bg r b{x) + —V 2 A(x) + ^V 2 b(x) 
ot Am Am 

-f\A(x)+g r b(x)\ 2 [A{x)+g r b(x)], 

ot \uf\ Am 

Here, the complex constants a — d and / can be obtained from the perturbation amplitudes of the path integral action. 
Their values depend on the propagators of noninteracting and free Fermi and Bose gases. These equations can be 
easily generalized to trapped gases as well. The dynamics represents a coupled system of Ginzburg-Landau equations 
(for the Fermi field A(x)) and Gross-Pitaevski equations (for the Bose field b(x)). If the atom molecule coupling g r is 
ignored, the dynamics of the Fermi and Bose fields decouple. This special case of the dynamics in Eq. [S] reproduces 
that of the single channel Fermi-BCS model, where the Cooper pairs do not form bound states, and the Feshbach 
detuning is small enough to keep them in the scattering continuum. In such a regime, Eqs. [5] decouple and the 
first equation becomes the time dependent Ginzburg Landau equation for the order parameter A(t). The nonlinear 
damping in the dynamics produces very rich and diverse behavior in Fermi-Bose mixtures. The dynamics of the 
single parameter Ginzburg-Landau equation, a special kind of Nonlinear Schrodinger equation, has been investigated 
thoroughly over several decades in various contexts ranging from Poiseuille flows (3l| , reaction-diffusion systems [32[ 
, as well as Fermi BCS superfiuids and BEC's [33[ . Stable solitonic solutions have been identified in certain regimes, 
as well as instabilities in lower dimensions obtained using various criteria ( for instance, Derrick's theorem [l^llHl or 
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the Vakhitov and Kolokolov stability criteria [35|, Hg| for U(l) invariant Hamiltonians) . These will not be reviewed 
here, and the reader is directed to ref [33| , and references therein, for details. 

General solutions to the dynamics in Eqs. [5] have been discussed by Machida and Koyama [2fj, as well as Chen 
and Guo [37j, Guo, Fang and Wang (38[, and others. The rest of this section focuses on the special case when 
quasiparticle excitations can be ignored. In that limit, the fermion field (characterized by the gap parameter) is 
spatially homogeneous due to it's coupling to the spatially homogeneous Bose field bo by momentum conservation. 
Then, the expansion of the effective action to the fourth order and subsequent computation of the functional derivatives 
with respect to the fields yields the mean field mean field dynamics of the order parameters following an analytic 
continuation of time to the real axis. Neglecting the excitations away from the BCS state of composite bosons in 
Eqs. [5] and projecting the dynamics into the space spanned by the BEC and fermion pair field yields the simplified 
dynamical system, 

*i + 47 (*i - * 2 ) ~ia^i+ z/3|*i| 2 *! = 
*2 + 2a*2 + 2ixl*2| 2 *2-iK7(*i-*2) = 0, (6) 

Here, ^j(t) = [A(t)Sis + 9rbo(t)] / (|^f|A/"), are the order parameter fields that describe the dynamics of the fermions 
in the BCS state and zero momentum bosons in the BEC respectively, and Af is the total (fermion) particle number. 
The constants expressed as Greek letters in equation [5] are evaluated from the perturbation amplitudes of the path 
integral. Their values depend on the propagators of a noninteracting Fermi gas in the trapping potential. These 
constants are complex and have been evaluated by Machida and Koyama [30| , and specific expressions for particle in 
a box confinement have been obtained by Roy [221 ] . 

The specific case of the particle in a box trap for the dual channel case, as covered by Roy (22j, is of particular 
interest, and can be easily generalized to other confinements. The equilibrium fixed points of the dynamical system 
in Eq. [6l together with the imposition of number conservation, J\f = lim | g_>. 00 <9 MF lnZ(/3) for the path integral grand 
partition function Z(j3), yields the BCS-BEC gap and number equations 

j(e F )(^ 1 -^ 2 )-a(e F )^ 1 +/3(e F )\^ 1 \ 2 ^ 1 = 0, 
2A,(e F )4' 2 +2x(eF)|*2| 2 ^2- K (e F ) 7 (e F ) (#i - # 2 ) = 0, 

2a 2 (e F )\y 2 \ 2 +e(tF)\*i\ 2 -l = 0. (7) 

Here, the new constant £ 2 = /j^d^a, e F = \i^f\, and the overbars denote the equilibrium values of the order parameters. 
These equations can be solved numerically to yield the initial conditions for the order parameters. Impulse quenches 
of the detuning v yield long term dynamics that can be approximated by linear displacements of ^i i2 about radially 
fixed trajectories obtained from Eqs. [7] viz. 

*i(t) = [r„ + J¥i(t)]e-**' t J 

* a (t) = [(l- J7 „)r I , + **i(t)]e-*"' t J (8) 



where r\ v — j^yfj^;, and r v = jj(a — 777^) is the equilibrium value obtained from Eqs. and the composite 



1/2 

bosons are presumed to be noninteracting (tig = 0) for the sake of simplicity. Substituting this into the dynamics 
in Eqns. |51 and retaining only terms up to linear order in ), yields a system of linear differential equations 

of the type 5ij ~ ^ fe Ji k 5xk- Here, 8^1,2 = &ci,3 + 16x2,4, and j,k = 1 — 4. The Jacobian matrix elements that 
describe this dynamics are given by ji k = dij/dxk at the fixed point. The linear dynamics is given by [IH, S3] 
\Sx(t)) ~ Ylj=i Cj- 1 f^^,) e°^- * , where Cj are constants, are the eigenvectors of J v that correspond to eigenvalues 
and \8x(t)) is the vector given by 8x\-^. The eigenvalues f2„ are given by the characteristic equation \ J V — = 0. 
In this case, the equation is biquadratic of the type ftf, -\-B u Vl 2 —C L , — 0, where the constants B v and C v are functions of 
a, (3, 7, k, and v [221 ] . The primary quantity of interest is the discriminant of the quartic equation 2?„ = 1 + (AC v /Bf^ . 

The roots of f2„ can be written in terms of this discriminant as VL V = ±y— B v /2 (l ± \fT>v) ■ The sign of T> v 
determines whether the roots are real or complex. Thus, it plays a key role in determining the nature of these 
trajectories. If f2„ are pure imaginary, a 'deep-BEC regime can be identified where the trajectories assume orbital 
motion about the radial fixed points. The orbits constitute two oscillations of frequencies f2 J = | yj— B v /2{1±^/V^) l / 2 \ . 
Their combination results in periodicity iff have rational winding numbers n + /n~ = f2+/r2~, where are integers 
with no common factors. The number density of the BEC, | <5\I> 2 1 2 = 8x\ + Sx 2 , will have 4 oscillations, 2f2^ and 
0+±f2~. When the fi„ start to pick up real parts, then a 'shallow BEC regime can be identified where the trajectories 
corresponding to negative real parts are stable (the eigenvalues with positive real parts correspond to displacements 
that violate number conservation and are unphysical) and decay in spirals to the orbits of the radial fixed points, thus 
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FIG. 1: (Color online) Plots of the signal |\I>2| 2 (i) from the solutions to Eqs. [6] for ub = after the system is quenched from a 
pure BEC to the shallow-BEC regime. Here, h = m = V = 1, and uf = —0.3. The left panels show plots of the time evolution 
for g r = 25.0 and v — —55.0, and the right panels show plots for g r = 40.0 and v — —140.0. The top panels show the time 
evolution for long times, and the bottom panels show the time evolution only after the onset of dynamics of C ( | <5\I/ 2 1 ) - The 
presence of partial collapse and revival of the matter wave can be noted after a relaxation time t r as indicated in the bottom 
panels. The revival times tn are also indicated in these figures. In the left panels, t r ~ 9, ta ~ 6, and in the right panels, 
t r — 4, tn ~ 2.5. Reproduced from [22| with permission of author. 



creating a stable limit cycle. Recently, Hopf bifurcations that are similar to those discussed above have been studied 
in the more general dynamics of Eqs.[S] [39(. 

A rapid quench from the 'deep BEC regime to the 'shallow BEC regime can thus be expected to yield orbital phase 
space dynamics of "J j around the equilibrium fixed points after a sufficiently long time, producing a 'collapse and 
revival' of the matter wave packet. The dynamics has been analyzed numerically for an ultracold gas of 100 fermions 
in a box of unit volume for larger values of in the shallow-BEC regime by Roy [22j and the ensuing dynamics of a 
representative parameter set are plotted in Fig. [TJ For long times, note from the panels in this figure that a sudden 
onset of partial collapse and revival of the matter wave begins after a short initial relaxation time t r . The onset of 
this collapse and revival takes place at around 10 3 numerical units. Also, note from these panels that all three time 
scales, td, t r and tn, are larger for smaller g r . Thus, it is expected that t& — > 00 as g r — > 0, removing the collapse and 
revival effect in the single channel case. 



III. MICROSCOPIC DYNAMICS 



This section details studies of the full microscopic dynamics of responses in a single channel fermionic superfluid in 
the BCS regime with time dependent parameters. The dynamics of the superconducting BCS state in metals has been 
of interest for a long time 41|, in particular, the regime where the BCS quasiparticles relax slowly compared to the gap 
(Anderson |42j). The regime of interest in this section is one where the external parameters can be changed rapidly 
compared to the equilibrium BCS relaxation time (HA~ 1 ) but slowly compared to the relaxation of the quasiparticle 
excitations. This allows for the analytic continuation in time of the equilibrium BCS problem. The resultant dynamical 
system consists of the Schrodinger equations of the microscopic wavefunction, and a time-dependent BCS gap equation 
that is self-consistently coupled to the Schrodinger equations. Among the time dependencies of interest are variation of 
the four fermion interaction amplitude via the Feshbach detuning, as well as variation of the chemical potential. Both 
of these are straightforward to perform in experiments. The four-fermion amplitude can be rendered time dependent 
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by varying the Feshbach field in time adiabatically (relative to the scattering continuum), and the chemical potential 
can be effectively varied in time by confining the system in a large isotropic trap with a time-dependent intensity, 
stiffness, or minima [13] . The Feshbach detuning of the confining lasers from the internal degrees of freedom of the 
confined atoms can also be varied in time to produce effectively time dependent chemical potentials [43[ 

In BCS systems modeled by self-consistent dynamics, both sudden and periodic quenches can be modeled by the 
BCS Hamiltonian 

H(t) = ^2[fk- fj,(t)] 4 ff c k <7 - g(t) ^2 c k+k»t c k'-k»i c k'i c kt- (9) 

kcr k,k',k" 

Here, /k = £k — Ho are the band energies relative to equilibrium chemical potential fig, and Cko- represent the 
annihilation operators for fermions of momentum k and spin a = {t?i}- The first term represents the kinetic energy 
of the fermions, and the second term the four-fermion interaction energy with amplitude g = g(0) > which represents 
attractive interaction between the fermions. In the rest of this work, the trap potential is assumed to vary sufficiently 
slowly that a locally constant chemical potential /io = €k F can be used to describe the fermions in the trap. At t = 0, 
this system can be mapped from the BCS Hamiltonian via a Hubbard- Stratonovich transformation (28l . [29| to yield 
a BCS-superfluid state from the Bogoliubov de-Gennes equations, 



«k \ _ ( (£k ~ Mo) A(k) \ f u k 
v u J { A*(k) -(e k -Mo) iU 



(10) 



where itk and are the amplitudes of the particle and the hole in a BdG quasiparticle and are related to the BCS 
wavefunction by 



nK+«kc kC Lk)io). (ii) 



The pair-potential A(k) depends on the pairing symmetry and is given by 



A(k) = Ao, s — wave, 

A(k) = Ao[cos(fca;) — cos(fcy)], d x 2_ y 2— wave. (12) 



For s-wave pairing, the pair potential satisfies the self-consistency relation 



A o = fY^WkVk- (13) 



k 



Eq. [TUl and IT3l admit the well-known BCS solution 



£(k) = ±V(£k-Mo) 2 + |A | 2 



W) = ^( 1 + (-)^)" 2 . (») 



As the system evolves in time, it can be modeled by the time dependent Bogoliubov de-Gennes equations viz. 

. R ( u k (t)\ _ ( [jk - fi{t)] A(k;t) \fuu(t) 
ldt {v k (t) ) ~ { A*(k;t) -[/ k -/*(*)] ){Mt) 

together with the self-consistency condition, which reads 



(15) 



A(k;*) ^ A(t) = g(t)Yu*Jt)v p (t), (16) 



p 



for s-wave pairing. Here, h = 1, and Uk and «k are the amplitudes of the particle and the hole in a BdG quasiparticle. 
Detailed justifications of the time dependent Bogoliubov de-Gennes equations can be obtained (in the reg ime of 
interest) from the microscopic nonequilibrium field theory of BCS systems. See, for instance, refs [161 ] and [17j | (Eqs. 
1-20)H| 
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The dynamics of the gap parameter in such regimes have been investigated by Szymanska et. al. [44J , Barankov et. 
al. (T3 |. Andreev et. al. [15|. Pekker ei. a/. (45j . using various trial solutions. The ensuing dynamics has, in particular, 
been cited as a way to identify pairing symmetries in Fermi superfiuids |24, 45} , as well as Stoner instabilities [46[ 
, and investigations of quench dynamics across the BCS-BEC crossover 46[ , where a renormalized g(t) switches 
sign for T = in the adiabatic limit 8]. The method used by Barankov, Levitov and Spivak (as well as a similar 
approach by Andreev, Gurarie and Radzihovsky [Hf) can be generalized in order to obtain the dynamics for fi(t) = 0, 
9(t) — 9o8(t) [H} ■ The unpaired state of no Cooper pairs, given by Uk(t) = 6*(/k)e~^ kt , t>k(i) = 6(— /k)e^ kt , satisfies 
Eqs. [15] with A(t) = Vt from Eq. [T5] Thus, this trajectory acts as a generalized fixed point of the dynamics in 
the form of a limit cycle. A linear stability analysis around this cycle was analyzed by Abrahams and Tsuneto [l3} 
, revealing that |A(i)| grows as e A °* for linear order displacements about the solutions of Uk(t) and t>k(t) above. 
Defining the variable % = (uk/2fk)(l + sig n /k) + (wk/2uk)(l — sign/k), the dynamics (for /k > 0) can be written as 

iiflk = 2 [/ k - /*(*)] w k + A(t) - A*(i)™ k , (17) 

with a similar equation for /k < 0. Substituting A(t) = e /a(t) into Eq. [17] and using the fact that Wk = 
2[/k - M(t)]/A* - id/dt(l/A*) yields 



2 + a* + |a| 2 (2ta£ - a) = 4 (a 2 - |a| 2 ) [/ k - m(*)] \a [/ k - mW] + «* ^ (18) 



The ansatz a(t) S R Vt, causes the RHS of Eq.[18] to vanish, yielding a nonlinear differential equation for the dynamics 
of a (and hence A) that is independent of k and decoupled from Eqs. [TS] , act — a 2 — I = 2ia 2 (i. In general, this leads 
to a contradiction, since this equation implies a(t) picks up an imaginary part, and relaxing the assumption a(t) G R 
causes momentum dependent terms in the RHS of Eq. [IS] to remain. Thus, this treatment does not yield a closed 
form for the dynamics of A(t). However, if //(f) = 0, Eq. [18]can be readily integrated with the ansatz a(t) £ R Vi [H} 
to yield A(t) = Aoe~' iSlt / cosh Aot. The frequency O can be obtained from the initial condition A(0) = Ao<?(0)/g 
(obtained by applying Eqs.[l4] to Eqs. [15] and the time derivative of Eq. [16] at t = 0) to yield ft = —ig(0)/g. Given 
that g{t) must be real for the BCS Hamiltonian to be Hermitian, this leads to the solution 



Ait) = 



A e- 



-r< 



cosh Ao£ ' 

(19) 



3, . . 



t=0 



If r > 0, the gap parameter starts from Ao and damps out to 0, stabilizing on the fixed point solution shown earlier 
in the para. However, if T < 0, then, for large times * > A^ 1 , A(i) - 2A e" (r+An)t -> if |T| < A , A(t) -> 2A if 
|T| = Ao, and A(t) blows up if |T| > Ao. The first of the three cases can be readily solved in the asymptotic limit. 
In the second case, Eqs. [TS] decouple from each other at large times with A(t) w 2Ao, and become non self-consistent 
linear Schrodinger equations with twice the equilibrium gap. The latter-most case is clearly unphysical due to norm 
conservation , |itk(f)| 2 + |fk(f)| 2 = 1 Vk and t. Investigations of this regime are a subject of further study. 

In the other two cases, the solution in Eq. [19] can be used to integrate Eqs. [17] . This yields w^it) — 
e -ia k (t) r^ k — i J dtA(t)e* ak (*)] for /k > 0, with a similar solution for /k < 0. Here, c?k is a constant determined 
from the initial conditions, and <2k(t) = J dt [/k — A{t)] can be evaluated in terms of Hyper-geometric polynomials. 
Analytical solutions to these equations have been provided by Barankov et. al. using Bloch vectors [l4| • The special 
case of a rapid quench in A deserves special interest. If the time that the gap takes to relax from Ao to the final value 
Ai (where Ai = if T > or — A < T < 0, and Ai = 2A if T = —A ) is neglected, then the final non self-consistent 

Schrodinger equations are those of the time independent Hamiltonian = TfEi(k), where Tk = T^x + r^y + t£z, 
with t£' v ' z being the three Pauli matrices in particle-hole space, and the vector Ei(k) = f^z + Aix with magnitude 
#i(k) = v/jfTAf. The initial state |^ k (0)) can be computed from the initial conditions for u k q and v k q given by 
Eqs. [H] . Thus, the solution to the particle-hole amplitudes in this limit is IV'k(i)) = exp (iH^t) |V>k(0)). This 
can be expanded using 5'C/(2) algebra, and can be used to determine the particle and hole amplitudes. If Ai = 0, 
then the system merely undergoes phase oscillations and the particle-hole amplitudes are frozen at their equilibrium 
values. If Ai = 2A , the particle-hole amplitudes can be averaged out over long times by averaging out the harmonic 
terms in their exact dynamics. Substituting the values of the equilibrium amplitudes from Eq. Q3] yields the long time 
average(s) to be |u k | 2 (|v k | 2 )i wnere 

h 





cK) = ^(l + H 7W=^ J ±TT^\ ^|1 + Ht7^I ■ (20) 
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FIG. 2: (Color Online) Numerical plots (blue and crossed) of Fermi surface magnetization after a double passage of the avoided 
crossing as a function of the drive frequency. The system parameters are Ao = fJ. a = 0.1, fio = 0.01, and the lattice size is 
N — 144 x 144 in two dimensions. This is contrasted with a plot of the analytical formula of irik F (red) obtained in section [Till 
Reproduced from [24[ with permission from the authors. 



Here, the approximation holds for large gap Ao. Thus, the system appears to attain a BCS-like steady state given 
by \ip s ) = rik( w k + w k c k c -k)|0); U P to a phase. If A is sufficiently large, then this state does not differ much 
from the equilibrium state as can be seen from the equation above. The steady state magnetization, defined by 
Q = 1 — 2 ^ k |t> k | 2 , is clearly about half of the equilibrium magnetization Q cq = 1 — 2 ^ k |t> k q | 2 for a sufficiently 
large gap. Thus, equilibrium states with half-filled symmetry retain this symmetry in the long time dynamics of the 
magnetization. 

The dynamics of impulse quenches in g can be formulated as a special case of the para above. Here, T — 0, and g is 
a constant for t > 0+. The system is populated in a BCS state characterized by a gap A = Ao. Barankov et. al. [14| 
analyzed the ensuing soliton dynamics for various values of Ao from Eq. [15] based on the linear stability analysis 
detailed above. This yields elliptic oscillations that cause overlapping solitons in U^Vk, causing them to oscillate 
about those for an ideal Fermi gas. Volkov and Kogan [17[ , have obtained generalized fixed points of [15] (with 
nonzero Ao), of the type it£uk ~ — *Ao0k, 1 ~ 2|«k| 2 ~ — i/k^k with real amplitude </>k, that restrict the Hilbert space 
to surfaces of constant u k i>k, i> k i>k- A linear stability analysis of Eqs. 1151 around this trajectory revealed bifurcations 
in the dynamics of A(t) that cause its trajectory to change from slow relaxations to oscillatory decays of rate 0(£ -1 / 2 ) 
towards Ao 7^ 0. The study of general time dependent g under these conditions is a subject of future study. 

The dynamics of and «k in the special case when g(t) is fixed at go — g, and fi(t) = fi a sin uit is the time periodic 
drive supplied to the system, has acquired renewed interest, and will be discussed below. In this regime, the analysis 
of Barankov et. al. detailed in previous paragraphs breaks down, necessitating other treatments for the ens uing 



dynamics. Such dynamics have been investigated in the context of Cooper-pair scattering by Challis et. al. 47 1 
using Bloch- Fourier series expansions for Uk,Vk and A(£), and truncating the series terms, obtaining approximate 
solutions that demonstrate Bragg scattering of correlated atom pairs. Studies of soliton formation in such systems 
at equilibrium due to the occurrence of Andreev reflection have also been performed by Antezza et. al. 48] , and 
the traveling dynamics of such solitons of the type A(x,t) — A(x — vt) for different pairing symmetries have been 
reported by Scott et. al. [4^| , Liao and Brand (5(| ■ Roy et. al. [24| have computed the response of the system as 
it evolves in time, demonstrating that the BCS self-consistency condition introduces nonlinearities in the dynamics 
and shapes the long-time behavior of the fermions subjected to the drive. If the system is started from the state 
in Eq. [14] at t — and evolved in time, a perturbative analysis of the ensuing dynamics is possible. The remaining 
part of this section focuses on the salient features of such a treatment done by Roy et. al. [24[ . The possibility of 
performing a linear stability analysis around a texture state similar to that by Volkov and Kogan as mentioned in the 
previous paragraph is a subject of future study. 

The primary response of interest is the effective 'magnetization' [231 ] . defined as m(t) = X^k m kWi wnere WkW = 
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FIG. 3: (Color online) Left Panel: Plot of Q, the time average of magnetization, as a function of u, with averaging carried 
over 10 drive cycles. Right panel: Plot of m(t) as a function of LJt/2-rr for representative values of cj indicated in the legend. 
A few representative values of the adiabatic magnetization m adb (t) is shown using crosses. In both panels, the equilibrium 
magnetization m(0) is indicated by a solid (red) horizontal line and all parameters are same as in Fig. [2] Reproduced from [24|] 
with permission from the authors. 
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FIG. 4: (Color online) Same as in Fig. [3] but for the non-self-consistent dynamics. Reproduced from [24 
the authors. 



with permission from 



( r k)W- The self-consistent dynamics given by Eqs. [T51 and ITol can be perturbatively integrated, using the Zener 
approximation j"> l] . for one passage across the avoided crossing in the adiabatic spectrum. The adiabatic state, is 
given by the equilibrium state in Eq. 1141 continued in time via the transformation fio — > no + fi a smujt, and the 
adiabatic spectrum by the energies ±E(k) transformed similarly. This integration has been done in a manner similar 

to [ll| except with a time dependent gap A(i). This yields a fermion density (which goes as |t>£ | 2 ) at the Fermi 
surface (given by /k = 0) after a single passage to be 



.,(1) |2 



l-2xo|c /k | 2 + 2y^xRe[(l + z)c /k 



(21) 



where 



f 7k 



oo 1 

T- 



(4 MaW ) Tl x d^t 



(22) 



t=o 



Here xo = ^0/(^)1 an d &w{t) = «k(i) x [A(^)/^o] x exp {i L dt'fi(t')}. For the non self-consistent case, A(t) = Aq 
Vi, and the series in Eq. |2"21 yields the traditional Landau Zener formula [5l|, |u 



(i) 12 



gap given by Eq. 1161 the quantities Ck and Iw^l 2 can be approximated by truncating Eqs 
Xo, yielding a fermion density that differs from that of the non self-consistent case. The density at the Fermi surface 



Xo /2. For a self-consistent 
^1and[2T1to lowest order in 
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after two complete passages (amounting to one period of the drive) across the avoided crossing can be approximated 
by [52[ l^ 2 '! 2 = 2|i£ 1 ' l | 2 ^1 — Iw^l 2 ^, where Stiickelberg interferometry has been neglected for large ui. The resultant 

density, written in terms of the 'effective magnetization' m^ F = (T^. =kp ) = 1 — 2|t£ 2 ' ) | 2 , approximately goes as xo/2- 
This yields a scaling law for this magnetization nik F ~ w -1 , that is similar to that of the non self-consistent case (as 
can be seen by Taylor-expanding the Landau Zener formula), but with a comparatively smaller scaling constant (by 
a factor of 4 to the lowest order). This result has been compared to numerical simulations in Fig. [5J The two diverge 
for to < Ao since that regime is nonperturbative, but agree quite well for larger lj. 

Numerical simulations of the dynamics above have been performed by Roy et. al. [24| for longer times while 
including the entire Brillouin zone in the calculation of the response. Plots of the instantaneous magnetization and 
its time average (over long times) in Fig [3] show many differences with the non self-consistent case, where Landau 
- Zener oscillations [52:] and other dynamical phenomena such as 'exotic freezing' of the response [2J| are expected. 
The primary difference is the appearance of a 'steady state' in the magnetization for the self consistent case when 
u) > Aq. Smaller w's yield results similar to the non self-consistent case, especially when u> -C Ao. In that regime, 
the magnetization approaches that of the adiabatic state ( which is defined in the fourth paragraph). The adiabatic 
magnetization, m adb , is shown as black crosses in the right panel of Fig. [3] This is consistent with the quantum 
adiabatic theorem, where adiabatic dynamics is expected if the time-scale of the dynamics (2n/u>) is greater than the 
relaxation time Ag -1 . 

The above behavior is to be contrasted with the non-self-consistent case summarized in Fig. |4j In the non self- 
consistent case, Eq. [13]is not analytically continued in time and the gap A(i) is kept constant at Ao. In this case, the 
dynamics is described solely by Eqs. 1151 without the self-consistent update in Eq.rTB]. Such dynamics can be described 
by the 1-D BCS reduced model [53| , whose co-ordinate space Hamiltonian 

H(t) = -£ [c\c l+1 + rcH +1 ] +^£(l- 24c) , (23) 

i i 

can be connected to the anisotropic quantum XY model. A Jordan Wigner transformation of this Hamiltonian at 
t = 54] , yields an SU(2) Hamiltonian similar to that in Eq. [10] with e& = — (/io/2) + cosfc and A(fc) = Tsinfc, 
where fj,Q = A*(0). The dynamics can thus be described by Eq. [T5] sans the self consistency in Eq. [T5] . This causes 
the system to decouple into Af independent two- level systems (henceforth abbreviated as TLS), each characterize d by 
momentum k. Such systems can also be generated from Ising (in 1-D, see Dziarmaga |5a|. Cherng and Levitov |56| . 
and Sen et. al, [STj) or Kitaev (in 2-D, see Chen and Nussinov [Hj], Sengupta et. al. [59{. Mondal et. al. [|3(|) models 
via the Jordan Wigner transformation, and their dynamics described by those of quantum TLS systems. 

Numerical solutions to such TLS dynamics for the isotropic case (A(k) = Ao) in 2-D are shown in Fig. [J] for the 
same parameters as those used for the self-consistent dynamics earlier in this section. Here, m(t) always executes a 
large, almost synchronized oscillation, approximately following the adiabatic path (black crosses in the right panel of 
Fig.rj}. Naturally, the resulting values of Q are close to m(0) (albeit with some small fluctuations). This suggests that 
the synchronous oscillation is simply a manifestation of the near-adiabatic nature of the dynamics. The qualitative 
departure from this behavior in the self-consistent case seems to stem form the non-adiabaticity induced by the self- 
consistency condition (Eq. II 6[) which makes the effective Hamiltonian non-linear. The overlapping eigenfunctions of 
the non-linear Hamiltonian makes the criteria for adiabatic behavior much more restricted compared to a linear case 
(see. e.g., Yukalov |61|). 

The quantum dynamics of driven TLS systems of the type described in Eq. [10] can be solved analytically using 
various techniques. The BCS solution in Eq. 1141 when continued in time via /xo — > ^{t), describe the solutions to 
Eq. [TU] in the adiabatic limit, when the characteristic time scale of the variation of fi(t) (w _1 if /i(t) = fi a s'muit) is 
much larger than the A -1 . In this limit, the system stays in the lower of the two adiabatic eigenstates E(k;t) — 
~ y/[fk — n(t)] 2 + A 2 (k). More rapid periodicities in a harmonically driven ji(t) may cause the system to vary across 
the quantum critical point (fiQ — 2), and produce defects away from the adiabatic state. The density of these 
defects can be computed for each k using Landau Zener Stiickelberg theory (reviewed in [52j ). In this limit, the 
dynamics of the amplitudes are presumed to be adiabatic for all times except when near an avoided crossing in the 
adiabatic eigenstates i.e when the separation between the two eigenstates is at a minimum. This occurs at times t c , 
given by the solution /k = fJ-(t c ) when the energy gap is 2|A(k)|. As the system varies from t = t c — to t — t c +, 
the nonadiabatic excitation density was estimated by Landau [62j and Zener [g3[ using the Zener approximation 
u(t) ~ /lawt around the avoided crossing. This process, illustrated by plots of E(k; t) in Fig. [5], yields a defect density 
of [5l| e _Xo /2, a result that is recovered from Eq. [21] in the non self-consistent limit. The Landau-Zener transitions 
contain classical signatures in the form of the Kibble-Zurek mechanism in the perturbative limit (see, for instance, 
Damski [64| , and Dziamarga [(35| ) . Even more rapid values of ui can produce dynamic symmetry breaking of the time 
averaged magnetization Q — limT^oo T~ x J Q dt m(t) as reported by Das [23[ . This occurs due to the extremely 
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FIG. 5: (Color online) Plots of the adiabatic energy levels ±E(k; i) as functions of dimensionless time illustrating the adiabatic 
impulse model in Landau-Zener theory in the non self-consistent case. The plots are for /k = 0, A(k) = Ao = /J, a /W, with 
fi(t) — fi a sinuit. The particle-hole amplitudes u^{t) and Vk(t) are expected to follow the adiabatic time dependence in Eqs. 1141 
with no continued in time to /io + fJ-(t), except at the avoided crossings (marked by arrows) when Landau-Zener transitions are 
expected. 
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FIG. 6: (Color online) Plot of the long time average of the magnetization Q vs r\ — 2/n a /w for multiple values of /x a for N = 100 
fermions subjected to the non self-consistent dynamics in Eg. 1151 with fj,(t) — /x a sinu;t, and ft = cosfc, A(fc) = sin A; for 
k G {— 7r,7r} (the Brillouin zone of a ID optical lattice with unit spacing). These are compared with the analytical result 
obtained using the Rotating Wave Approximation in [2^|. The peaks Pi, P2 and P3 represent maximal freezing, corresponding 
to three zeros of Jo(??), occurring at r\ — 5.520..., 8.653..., and 11.971... respectively. Reproduced from [23|] with permission from 
author. 



nonadiabatic nature of the dynamics if u> 3> A. This also produces a dynamical hysteresis similar to that from the 
Landau Zener dynamics. Das [23| derives the onset of such hysteresis using the Rotating Wave Approximation ( 
henceforth abbreviated as RWA) for large lj in a manner similar to Oliver et. ah [6q . as well as Ashhab et. al. [67| 
(also see 52] ) . Eq. [15] can be approximated by an RWA Hamiltonian that contains only the slowest - oscillating 
contribution to the time translation operator, given by Uk(t) = exp {ir z J dt [/k — ^(t)]} in the interaction picture. 
This can be obtained by a Fourier series decomposition of a periodic Uk(t) and ignoring all the oscillating terms (23j . 
The resultant time independent RWA Hamiltonian is simply #™ a ~ A(k) Jo(r])T v , where 77 = 2fi a /uj for amplitude fi a , 
and Jo(x) is the zeroth order Bessel function. The Schrodinger equation for this Hamiltonian can readily be integrated 
to yield solutions where the amplitudes itk and start to show hysteresis for parameters that approach Jo (77) = 0, 
and dynamically freeze at their equilibrium values at the zeroes of Jo (77), i.e when 7J rwa vanishes. This can be seen 
in the numerical simulations plotted in Fig. [5] for 100 fermions in a 1-D optical lattice of unit spacing. Here, the 
pairing symmetry is anisotropic, and only the non self-consistent solutions to Eq. [TS] are evaluated. The comparison 
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with the results obtained from the RWA are excellent even at low w, but are expected to differ when uj < A, where 
either the adiabatic or Landau-Zener formalisms yield better results. Unlike the Landau-Zener limit, however, the 
dynamical hysteresis follows a mechanism that is very different from the quantum Kibble- Zurek mechanism [23| and 
thus illustrates a uniquely quantum phenomenon. 

IV. CONCLUSIONS 

This article, has reviewed numerous works investigating the dynamics of quantum quenching in two distinct regimes 
of Fermi BCS systems, classified according to the time scales of the quench relative to the equilibrium BCS gap. In one 
case, the system can be described in the mean field by the time dependent Ginzburg Landau equations of the fermion 
gap parameter, with a coupling to the Gross Pitaevski Bogoliubov equations for molecular Bosons. Such systems 
admit to stable solitons, as well as unstable ones in lower dimensions, and impulse quenches away from equilibrium 
can cause collective behavior that manifest as collapses and revivals of the matter waves. In the other case, the 
system can be described by self-consistent Bogoliubov de-Gennes equations whose dynamics produce behavior that 
is markedly different from the non self-consistent dynamics that is appropriate for quenching in Ising and Kitaev 
models. 

The author thanks CSIR, India for support under the Scientists' Pool Scheme No. 13(8531— A)/2011/Pool. Grat- 
itude is also extended to Prof. K. Sengupta, IACS Kolkata, as well as Dr. A. Das, Max Planck Dresden for their 
critiques. The author also thanks Prof. G Ambika, USER Pune and convener of the NCNSD conference, 2012, for 
inviting this review, as well as Prof. A. Laxminarayan, IIT Madras, for convening the quantum chaos minisymposium 
at the NCNSD conference, the proceedings of which form a part of this review. 
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